/*
 * Class:        InverseGaussianProcess
 * Description:  
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       
 * @since

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */
package umontreal.ssj.stochprocess;
import umontreal.ssj.rng.*;
import umontreal.ssj.probdist.*;
import umontreal.ssj.randvar.*;

/**
 * The inverse Gaussian process is a non-decreasing process where the
 * increments are additive and are given by the inverse gaussian
 * distribution,  @ref umontreal.ssj.probdist.InverseGaussianDist. With
 * parameters @f$\delta@f$ and @f$\gamma@f$, the time increments are given
 * by  @ref umontreal.ssj.probdist.InverseGaussianDist @f$(\delta dt/\gamma,
 * \delta^2 dt^2)@f$.
 *
 * [We here use the inverse gaussian distribution parametrized with
 * IGDist@f$(\mu,\lambda)@f$, where @f$\mu=\delta/\gamma@f$ and
 * @f$\lambda=\delta^2@f$. If we instead used the parametrization
 * @f$IGDist^{\star}(\delta, \gamma)@f$, then the increment distribution
 * of our process would have been written more simply as
 * @f$IGDist^{\star}(\delta dt, \gamma)@f$.]
 *
 * The increments are generated by using the inversion of the cumulative
 * distribution function. It therefore uses only one
 * @ref umontreal.ssj.rng.RandomStream. Subclasses of this class use
 * different generating methods and some need two
 * @ref umontreal.ssj.rng.RandomStream ’s.
 *
 * The initial value of this process is the initial observation time.
 *
 * <div class="SSJ-bigskip"></div><div class="SSJ-bigskip"></div>
 */
public class InverseGaussianProcess extends StochasticProcess {

    protected RandomStream       stream;

    protected double   delta;
    protected double   gamma;

    protected double   deltaOverGamma;
    protected double   deltaSquare;
    // mu and lambda are the common names of the params for InverseGaussianGen.
    protected double[] imu;
    protected double[] ilam;

    // Number of random streams needed by the current class
    // to generate an IG process.  For this class, = 1, for subclasses
    // will sometimes be, = 2.
    int numberOfRandomStreams;

    // needed by InverseGaussianProcessMSH
   protected InverseGaussianProcess()  { }

/**
 * Constructs a new `InverseGaussianProcess`. The initial value @f$s0@f$ will
 * be overridden by @f$t[0]@f$ when the observation times are set.
 */
public InverseGaussianProcess (double s0, double delta, double gamma,
                                  RandomStream stream) {
        this.x0 = s0;
        setParams(delta, gamma);
        this.stream = stream;
        numberOfRandomStreams = 1;
    }


   public double[] generatePath() {
        double s = x0;
        for (int i = 0; i < d; i++) {
            s += InverseGaussianGen.nextDouble(stream, imu[i], ilam[i]);
            path[i+1] = s;
        }
        observationIndex   = d;
        observationCounter = d;
        return path;
    }

/**
 * Instead of using the internal stream to generate the path, uses an array
 * of uniforms @f$U[0,1)@f$. The array should be of the length of the number
 * of periods in the observation times. This method is useful for
 * @ref NormalInverseGaussianProcess.
 */
public double[] generatePath (double[] uniforms01) {
        double s = x0;
        for (int i = 0; i < d; i++) {
            s += InverseGaussianDist.inverseF(imu[i], ilam[i], uniforms01[i]);
            path[i+1] = s;
        }
        observationIndex   = d;
        observationCounter = d;
        return path;
    }

   /**
    * This method does not work for this class, but will be useful for the
    * subclasses that require two streams.
    */
   public double[] generatePath (double[] uniforms01, double[] uniforms01b) {
       throw new UnsupportedOperationException("Use generatePath with 1 stream.");
    }


   public double nextObservation () {
        double s = path[observationIndex];
        s += InverseGaussianGen.nextDouble(stream,
                  imu[observationIndex], ilam[observationIndex]);
        observationIndex++;
        observationCounter = observationIndex;
        path[observationIndex] = s;
        return s;
    }

/**
 * Sets the parameters.
 */
public void setParams (double delta, double gamma) {
        this.delta = delta;
        this.gamma = gamma;
        deltaOverGamma = delta/gamma;
        deltaSquare    = delta*delta;
        init();
    }

   /**
    * Returns @f$\delta@f$.
    */
   public double getDelta() {
        return delta;
    }

   /**
    * Returns @f$\gamma@f$.
    */
   public double getGamma() {
        return gamma;
    }

   /**
    * Returns the analytic average which is @f$\delta t/ \gamma@f$, with
    * @f$t=@f$ `time`.
    */
   public double getAnalyticAverage (double time) {
        return delta*time/gamma;
    }

   /**
    * Returns the analytic variance which is @f$(\delta t)^2@f$, with
    * @f$t=@f$ `time`.
    */
   public double getAnalyticVariance (double time) {
        return delta*delta*time*time;
    }


    protected void init () {
        super.init(); // set path[0] to s0
        if(observationTimesSet){
            x0 = t[0];
            path[0] = x0;
            imu  = new double[d];
            ilam = new double[d];
            for (int j = 0; j < d; j++)
            {
                double temp = delta * (t[j+1] - t[j]);
                imu[j]  = temp / gamma;
                ilam[j] = temp * temp;
            }
        }
    }

   public RandomStream getStream () {
        // It is assumed that stream is always the same
        // as the stream in the InverseGaussianGen.
        return stream;
    }

   public void setStream (RandomStream stream) {
       this.stream = stream;
    }

/**
 * Returns the number of random streams of this process. It is useful because
 * some subclasses use different number of streams. It returns 1 for
 * @ref InverseGaussianProcess.
 */
public int getNumberOfRandomStreams() {
        return numberOfRandomStreams;
    }

}